#include <iostream>
#include <iomanip>
#include <sstream>
#include <cstdlib>
#include <map>
#include <algorithm>
#include <fstream>
#include <time.h>

#include "fastjet/ClusterSequence.hh"
#include "fastjet/tools/Pruner.hh"

#include "Njettiness.hh"

#include "CmdLine.hh"
#include "TFile.h"
#include "TChain.h"

#include "TH1F.h"
#include "TH2F.h"

#include "MyAnalysis.h"
#include "subAnalysis.h"

int main(int argc, char** argv)
{

  //initialize command line interface
  CmdLine cmdline(argc,argv);

  if(!cmdline.present("-i") || !cmdline.present("-o"))
    {
      cout<<"Usage: -i <input> -o <output> -a [algorithm] -r [radius] -m [mass] -n [nev] -p -w -gluino -rutgers"<<endl;
      cout<<"-a \t clustering algorithm to be used"<<endl;
      cout<<"-r \t clustering radius"<<endl;
      cout<<"-m \t a mass scale for rescaling the plots"<<endl;
      cout<<"-n \t number of events to process"<<endl;
      cout<<"-p \t whether jet pruning is used"<<endl;
      cout<<"-w \t whether event has a weight stored as evt Number"<<endl;
      cout<<"-gluino \t retrieve plots that trace the gluinos "<<endl;
      cout<<"-jet4 \t apply 4 jets cut "<<endl;
      cout<<"-jet2 \t apply 2 jets cut "<<endl;
      cout<<"-display \t number of display plots"<<endl;
      cout<<"-rutgers \t enable rutgers analysis plots"<<endl;
      

      return 1;
    }

  //create analysis object
  MyAnalysis myanalysis(cmdline);
  myanalysis.tree_init();
  myanalysis.hist_init();
  myanalysis.hist_init2();
  myanalysis.analyze();

  return 0;
}


//main function to loop over all the entries
void MyAnalysis::analyze()
{
  int n_acc_event=0;
  //begin event loop
  for(int entry=0; 
      (maxevt<0 || n_acc_event<maxevt) && entry<tree.GetEntries();
      entry++)
    {

      if(entry%1000==0 && entry/1000 > 0)
	cout<<"MESSAGE: "<<entry<<" Processed"<<endl;

      //populate our variables
      tree.GetEntry(entry);

      Event_=entry;
      if(cmdline.present("-unit-weight"))
	{
	  weight=1.0;
	}

      if(weight==0.0)
	weight=1.0;

      //grab all the jets
      makejets();
      
      if(entry < 10)
	{
	  cout<<"Event "<<entry+1<<": "<<endl;
	  output();
	}

      //fill_rmass();

      if(!cut()) continue;
          
      n_acc_event++;

      //fill trees for "event display" plots    
      fill_display();

      //fill masspt;
      fill_masspt();
      
      //fill Nsubjettiness
      fill_Njettiness();
            

      /*
      if(cmdline.present("-jet4"))
	{
	  fill_rutgers();
	}
      */


    }

  mass_profile->Write();
  fileout->Write();
  fileout->Close();
}



void MyAnalysis::analyze_glu(int i_key)
{
  if(!b_gluino)
    return;
  using namespace fastjet;

  //grab indices for gluinos
  int iglu1=0, iglu2=0;
  for(int i=0; i<size; i++)
    if(PID[i]==1000021 && Status[i]==3)
      {
	if(iglu1==0)
	  iglu1=i;
	else
	  {
	    iglu2=i;
	    break;
	  }
      }

  //grab the first status 3 decay products
  set<int> glu1_index1, glu2_index1;
  for(int i=0; i<size; i++)
    {
      if(M1[i]==iglu1)
	glu1_index1.insert(i);
      
      else if(M1[i]==iglu2)
	glu2_index1.insert(i);
          
      if(Status[i]==2)
	break;
    }

  //  assert(glu1_index1.size()==3);
  //  assert(glu2_index1.size()==3);
  

  //grab the first status 3 decay products
  set<int> glu1_index2, glu2_index2;
  for(int i=0; i<size; i++)
    {
      if(glu1_index1.count(M1[i]))
	glu1_index2.insert(i);
      
      else if(glu2_index1.count(M1[i]))
	glu2_index2.insert(i);         
    }

  PseudoJet g1,g2;
  //grab the first status 2 decay products
  for(int i=0; i<size; i++)
    if(Status[i]==2 && PID[i]!=92)
    {
      if(glu1_index2.count(i))
	{
	  PseudoJet j(Px[i], Py[i],
		      Pz[i], E[i]);
	  hist2[key_gluino1[i_key]]->Fill(j.eta(),
					  j.phi(),
					  j.pt());
	  
	  g1+=j;
	}
	
      else if(glu2_index2.count(i))
	{
	  PseudoJet j(Px[i], Py[i],
		      Pz[i], E[i]);
	  hist2[key_gluino2[i_key]]->Fill(j.eta(),
					  j.phi(),
					  j.pt());
	  
	  g2+=j;	  
	}

          
    }
 
}


